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We present large statistics simulations of 3-dimensional star polymers with up to / = 80 arms, 
and with up to 4000 monomers per arm for small values of /. They were done for the Domb- Joyce 
model on the simple cubic lattice. This is a model with soft core exclusion which allows multiple 
occupancy of sites but punishes each same-site pair of monomers with a Boltzmann factor v < 1. 
We use this to allow all arms to be attached at the central site, and we use the 'magic' value v = 0.6 
to minimize corrections to scaling. The simulations are made with a very efficient chain growth 
algorithm with resampling, PERM, modified to allow simultaneous growth of all arms. This allows 
us to measure not only the swelling (as observed from the center-to-end distances), but also the 
partition sum. The latter gives very precise estimates of the critical exponents 7/. For completeness 
we made also extensive simulations of linear (unbranched) polymers which give the best estimates 
for the exponent 7. 



I. INTRODUCTION 

Star polymers are of interest both for their technical 
applications, ranging from lubricant additives to paints 
0, 13 , and for the theoretical challenge which they repre- 
sent. Polymer theory in general is one of the prime fields 
where renormalization group theory can be used and 
compared in detail with real experiments 3, 4J. The sim- 
plest non-trivial objects in this respect are the partition 
sum and the r.m.s. end-to-end distance of a single long 
flexible linear (unbranched) polymer with A" monomers 
in a good solvent, which scale as 
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Star polymers, i.e. / such chains linked together at a sin- 
gle point, are some of the simplest examples of polymers 
with non-trivial topology. As shown by Duplantier 
all such polymer networks are characterized by equations 
similar to Eqs. and J2J, with the critical fugacity \x 
and the critical exponent v being the same for all topolo- 
gies, but with 7 being universal only within each topol- 
ogy. For star polymers composed of / arms of length A" 
each, one has in particular 
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where Rnj is the r.m.s. Euclidean center-to-end dis- 
tance. 

The behaviour of 7/ and of the swelling factor Af are of 
central interest, both for finite / and for / — ► 00. In two 
dimensionSj_7/ can be calculated exactly using conformal 
invariance pg, but no exact results are known for d = 3. 
Renormalization group methods give e expansions up to 



third power in e = 4 — d jg , but these are non-convergent 
power series and have to be resummed before being appli- 
cable in d — 3. The results are debated, in particular for 
large values of / |7j • For the swelling factor the situation 
is similarly unclear. Phenomenologists tend to compare 
with predictions based on Gaussian (i.e. free) chains @ 
or on heuristic assumptions 0, El • There exist several 
renormalization group calculations, but those not based 
heavily on simulation data 0,0] seem to describe some 
of the data rather poorly, and Monte Carlo simulations 
are needed to fix free parameters in such theories E|Jp1 • 

In view of this, Monte CarloEl E El [H EMlla 
and molecular dynamics [2ll |22| simulations have played 
a major role in the efforts to understand the behaviour of 
star polymers. Molecular dynamics simulations |22j| have 
indeed been used to study very large stars, with up to 80 
arms of length A" = 100 each, but it is not clear whether 
these simulations have really reached equilibrium. More- 
over, both molecular dynamics and Monte Carlo methods 
with fixed chain lengths (including the pivot algorithm 
0>H3) cannot measure the partition sum and thus give 
no information on . For the latter one has to use chain 
growth methods Il5l 13 IH EHH E3. Unfortunately, 
with the methods used so far it has not been possible to 
go beyond 24 arms [l?]], and even these were too short 
and the data were too noisy to provide a clear cut picture 
of the asymptotic behaviour. 

We decided therefore to perform simulations with sev- 
eral improvements which allow us to reach much larger 
systems and much higher accuracy. To obtain a good es- 
timate for fi and for the critical exponents of unbranched 
polymers, we also made extensive simulations of linear 
chains. The model and the method of simulation are de- 
scribed in the next section. Results are given in Sec. 3, 
while we end with a discussion in Sec. 4. 



II. MODEL AND METHOD 

Let us first describe in detail our model. For efficiency, 
and since we are only interested in scaling behaviour, we 
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use a lattice model. Indeed, we use the simplest ver- 
sion, the simple cubic lattice. But instead of simulating 
self avoiding walks as in previous works, we simulate the 
Domb- Joyce model [24| at its 'magic' interaction strength 
v — v* . In the Domb- Joyce model polymers are described 
by lattice walks where monomers sit at sites and are con- 
nected by bonds of length 1 . Multiple visits to the same 
site are allowed, but for any pair of monomers occupying 
the same site one has a repulsive energy e > giving rise 
to a Boltzmann factor v — exp(— j3e) < 1. The partition 
sum of a linear chain molecule of N + 1 monomers is thus 
a sum over all walks of N steps, each weighted with v m 
where m is the total number of pairs occupying the same 
site, m = ^i < j ^x^Xj- For star polymers we studied in 
the present work two variants. In both variants arms of 
N monomers are attached to a central site. In the first 
variant, the central site is singly occupied. In the sec- 
ond, it is occupied by / monomers, one for each arm. 
We studied both variants in order to verify that results 
were independent of this detail, and we include in our 
final error estimates the uncertainty it entails. 

Using the Domb- Joyce model has two main advan- 
tages. First of all, it allows us to attach a large number 
of arms to a point-like center. In the present work, we go 
up to / = 80 [2j| . Previously, authors had used extended 
cores. Although these cores were much smaller than the 
radii of the polymers themselves and should thus not de- 
stroy the asymptotic scaling, they do introduce a finite 
length scale and present therefore corrections to scaling 
terms which complicate the analysis. 

More important is that there is one special ('magic') 
value of v, called v* in the following, where corrections to 
scaling are minimal and where asymptotic scaling laws 
are reached fastest. For single chains it has been esti- 
mated p6l |27| as v* « 0.6 with rather small error bars, 
and we shall in the following assume this value to be ex- 
act. In the renormalization group language, the flow of 
the effective Hamiltonian to its fixed point in the stable 
manifold of the latter contains one direction of slowest 
approach. For a generic starting point there is a non-zero 
component in this direction, which then determines the 
leading correction to scaling. If one starts however with 
the flow such that this component is absent, the approach 
to scaling is governed by the next-to-leading correction 
term and is much faster. A similar observation has been 
made also for off-lattice bead-rod models with fixed bond 
length, where the leading corrections to scaling are ab- 
sent for a certain 'magic' ratio between bead size and rod 
length [2l|2i|I3. 

Since the value of v* should depend only on the inter- 
nal structure of the chains, for star polymers it should 
be independent of the number of arms. For the bead- 
rod model this was carefully verified in 14]. We thus 
simulated only with v — v*. 

For our simulations we used the pruned-enriched 
Rosenbluth method (PERM) 0] . This is a biased chain 
gro wth algorithm, similar to the Rosenbluth-Rosenbluth 
[32| method. In the latter, the bias induced by avoid- 



ing double occupancy is compensated by a weight fac- 
tor which is basically of entropic origin. In the present 
case, we have both a bias compensating factor and a 
Boltzmann factor, the product of which tends to fluc- 
tuate wildly if there is no perfect importance sampling. 
These fluctuations are suppressed in PERM by pruning 
low weight configurations and cloning those with high 
weight. Indeed, any bias can be employed in PERM, as 
long as it can be compensated by a weight factor. In pre- 
vious simulations of diluted polymers we use a Markov 
approximation called Markovian anticipation [33l. l34l.l35j . 
In the present case we did not expect this to be very use- 
ful, because the main interactions are not within one arm 
but between different arms. Thus we used instead a very 
simple bias where each arm tends to grow preferentially 
outward (except for the simulations for / = 1 where we 
used of course Markovian anticipation). The strength of 
this bias was adjusted by trial and error. It decreased 
with the length of the arm and increased with /. De- 
tails will not be given since they are not very important, 
and working without this bias would have increased the 
errors by only a factor « 2, in general. 

A final comment is that it is easy to modify the basic 
PERM algorithm given e.g. in the appendix of [3l| such 
that all / arms are grown simultaneously |36j . This is 
done by having / growth sites Xi, . . . x/. Chain growth 
is made in PERM by calling recursively a subroutine for 
each monomer addition. For multi-arm growth, we add 
an integer k S [1, . . . /] to the argument list of this sub- 
routine, such that a subroutine called itself with argu- 
ment k calls the next subroutine with (k mod /) + 1. 
In this way a monomer is added to each arm before the 
next round of monomers is added. Compared to a scheme 
where one arm is grown entirely before the next arm is 
started, the main advantage is that each chain grows in 
the field of all the others, and is thus, by the population 
control (pruning/cloning), guided to grow into the cor- 
rect outward direction. If chains were grown one after 
the other, this bias would be absent for the first chains 
which then would grow into "wrong" directions, resulting 
in very low weight configurations. 



III. RESULTS 
A. Partition Sums and 7-Exponents 

One of the outstanding features of chain growth meth- 
ods such as PERM is that they give estimates for the 
partition sum. Indeed, these estimates are a basic part 
of the simulations, since the population control is based 
on these estimates. 

According to Eq. (3) we expect Znjh* n to approach 
a power law const A™ -1 at large N. One precise 
way to estimate 7/ is to subtract a term afhiN from 
ln(.£jv,//-t ), an d adjust the constant a/ such that the 
difference gives a flat curve for large N, when plotted 
against InN. This gives then 7/ = 1 + a/. Alternatively, 
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we could plot In Znj — In Zjn,i — off In N against In AT, in 
which case a flat curve is obtained when a' — 7/ — 7. We 
prefer both methods over a least square fit, say, since they 
allow directly to check visually for the presence of correc- 
tions to scaling. If such corrections seem needed, one can 
subtract them and obtain in this way the most reliable 
and precise estimates of 7/ - remembering of course that 
estimating a critical index involves an extrapolation and 
is thus ill-posed, giving at best subjective error estimates. 

For either method we need precise estimates of the par- 
tition sum of linear chains. We thus performed first ex- 
tensive simulations of linear (/ = 1) Domb- Joyce chains, 
creating altogether w 4 x 10 s chains of length AT = 8000. 
In Fig. 1 we plot effective exponents obtained from triple 
ratios (2|| Z^ n Z^ n /Zn- Here a and b are chosen such 
as to minimize statistical and systematic errors |26). and 
powers x and y are fixed such that fj, and the overall 
normalization drop out. With a — 1/3 and b — 5 we 
have 



71nZ/v — 61nZ 



JV/3 



InZ, 



ln(3 6 /5) 



(5) 



which is plotted against 1/N - 96 . The fact that we find 
essentially a straight line (apart from odd/even oscilla- 
tions due to the special structure of the cubic lattice) 
indicates that the leading correction to scaling exponent 
is A ?s 0.96 which is much larger than the value Awl/2 
for generic self avoiding walks, indicating that v = 0.6 is 
indeed close to the magic value. Our estimate is therefore 
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This is in good agreement with the best previous esti- 
mates [26l [37J but more precise. Using this estimate we 
obtain then 



(i = 0.18812145 ± 0.00000003. 
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After having obtained a precise estimate for fj,, we can 
now discuss the results for stars. Results for a few se- 
lected values of / are shown in Fig. 2. We plot there 
In Znj + /AHn/z for both variants, i.e. the center singly 
occupied or / times occupied. The latter gives smaller 
values of Znj, but the difference is visible only for 
N = 2. For larger N both agree, except for / = 80 and 
large values of N where our sampling algorithm starts to 
break down. 

For a precise estimate of 7/ we of course did not use 
plots like Fig. 2, but we subtracted (7/ — 1) In AT as ex- 
plained above. Then we see (Fig. 3) that there are non- 
negligible corrections to scaling, but our arms are long 
enough so that our estimates of 7/ are not affected by 
them. Our final results, obtained by averaging over both 
variants of the model, are shown in Table 1 and in Fig. 4. 
In Table 1 we also give additional information such as the 
arm lengths and the total statistics. We also list previous 
estimates for comparison. We see reasonable agreement 
in general, although those previous estimates which were 
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FIG. 1: Effective exponents y e g(N) for linear (/ = 1) 'magi- 
cal' Domb- Joyce polymers against 1/N 0,96 . 
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FIG. 2: Logarithms of the partition functions Znj multiplied 
by n* . For each pair of close- by lines, the upper refers to 
singly occupied centers, while the lower one has / monomers 
located at the center. 



quoted with error bars pjj are off by many standard 
deviations. We should add that the simulations in [Tr| 
involved much shorter chains and lower statistics. 

Previous theoretical predictions of 7/ used e = 4 — d 
- expansions [|| Q and the cone approximation [3^, |3t| . 
The latter assumes that each branch is confined to a cone 
of space angle 47r//, and gives 
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As seen from Fig. 4 this is not too far off, but it definitely 
does not provide a quantitative fit to our data. The best 
fit with a power law 7/ — 1 ~ — (/ — 1.5) 2 would be 
obtained with z 1.68, but we do not claim that this 
exponent has any deeper significance. 

In contrast to the cone approximation which is basi- 
cally heuristic and cannot be improved systematically, e - 
expansions have a firm theoretical basis. But the expan- 
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TABLE I: Main results. 
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FIG. 3: Logarithms of the partition functions Zj?j multiplied 
by fi^ N , plus (1 — 7/) In N. For each pair of close-by lines, the 
upper refers again to singly occupied centers, while the lower 
one has / monomers located at the center. 
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FIG. 4: Exponents 7/ versus /. The full line is just a polygon 
connecting the points, the dashed line is a fit with the large-/ 
behaviour as predicted by the cone approximation, Eq. JHJ . 
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FIG. 5: Effective exponents v e ff(N) for linear 'magical' Domb- 
Joyce polymers against 1/iV ' 96 . 



sion itself is at best asymptotic, and each term gives a 
contribution to jf which is polynomial in /. Thus it can- 
not be used without re-summation. Such re-summations 
have not yet been attempted for / — ► oo. For small /, 
results are given in 0,0, and are listed in Table 1. They 
are in the right order of magnitude, but their precision 
is not sufficient to draw any firm conclusion beyond the 
fact that the resummed e - expansion is obviously not in 
conflict with the Monte Carlo data. 



B. Coil Sizes 

We measured only rms. center-to-end distances of the 
arms (resp. end-to-end distances for / = 1). This was 
done 'on the fly', i.e. we did not store each configuration 
and measure its properties in a second step of analysis. 
The reason is that an off-line analysis would have re- 
quired very large files, and reading a configuration from 
disk or tape would have been not much faster than cre- 
ating it from scratch. We neither measured shape pa- 
rameters nor radii of gyration, since any such additional 
measurement would have slowed down the analysis con- 
siderably, and since the main purpose of the present work 
was to demonstrate the efficiency of PERM and to study 
the main universal properties of large stars. 

As for jf, we first need a careful analysis of lin- 
ear chains, to obtain precise estimates of v and of the 
amplitude A\. In Fig. 5 we plot effective exponents 
v eS (N) = (h3.16)- 1 hx[H^ N /R 2 N ], again versus 1/iV - 96 . 
We see again a straight line (as in Fig. 1), verifying again 
that the correction to scaling exponent is close to 1. Ex- 
trapolating to N -> oo we find v = 0.58767(20). To- 
gether with previous estimates reviewed in |35| , this gives 
our best estimate 

v = 0.58765 ±0.00020 . (9) 

Notice that this is more precise than the field theoretic 
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FIG. 6: R 2 NJ /N 2v plotted versus In N, for seven selected val- 
ues of /. For each /, results are shown for both variants of 
the model (single occupancy of the center: lower curve; /-fold 
occupancy: upper curve). The structures seen for / > 40 and 
large N are statistical fluctuations. 



estimates obtained from the e-expansion Q. The result- 
ing amplitude is then 

Ai = lim R 2 N /N 2v = 0.8038 ± 0.0005 . (10) 

N^oo 

To obtain the amplitudes Af for stars we assume 
the value of v as given above. We can then plot ei- 
ther R 2 Nf /N 2v versus N (which gives Af directly), or 
R 2 N j/R% versus N, which gives Af/A\. To check for 
systematic corrections we did both. Some typical curves 
obtained with the first method are shown in Fig. 6. For 
each value of / we see two curves, one for each variant 
of the model: The upper curve is always that with the 
center / times occupied, the lower one corresponds to a 
singly occupied center. For large values of / (/ > 40) we 
see large fluctuations, indicating the limit where our sam- 
pling breaks down. Otherwise we see large corrections to 
scaling, but they all are dominantly ~ 1/N, i.e. analytic 
corrections, and they have rather small influences on our 
final estimates of Af. 

These estimates are given in Table 1 and plotted in 
Fig. 7. We show indeed the ratios Af/A\, to facilitate the 
comparison with previous estimates. The best previous 
estimates are those of Zifferer [2(J and are also given in 
Table 1. We see very good agreement, even if most of the 
values of [2(j are outside our error bars. The data of Zif- 
ferer were obtained from simulations on the tetrahedral 
lattice, and they indicate that the ratios A t jA\ are in- 
deed universal. In , stars with up to 80 arms were sim- 
ulated off-lattice by means of molecular dynamics. But 
it seems that the stars with N — 100 were not equili- 
brated, at least for / = 1 (see Table 1 of l^)- Therefore 
we list for comparison only the data for N = 50. They 
are systematically larger than our results and those of 
[2fj| . and the discrepancy increases with /. This suggests 
that even these simulations had not reached equilibrium. 
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FIG. 7: Amplitude ratios Af/Ai plotted against /. 



The most cited predictions for At I Ax are from a 
heuristic blob model 0, . It gives 
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which is in gross violation with our data. The fact that 
this Daoud-Cotton model gives a too strong swelling with 
/ is well known [23, E3 ■ Our data cannot be fitted by 
a pure power law, but asymptotically (for / — > oo) they 
tend roughly to AfjA\ w 0.78 /°- 235 . Again we do not 
expect this to be the true asymptotic behaviour, but it 
provides at least a useful guide for extrapolations. 

Renormalization group (RG) calculations of star poly- 
mer sizes have been performed in 0, 0, but it was 
already pointed out in 0, 0] that these have difficul- 
ties in describing the large-/ behaviour. Using their own 
simulations to fix some of the parameters in an improved 
RG calculation, Lue and Kiselev 0, 0] were able to 
fix these problems in the sense that their RG calcula- 
tion described perfectly the behaviour of the penetration 
function |l3ul4j . Unfortunately, they did not give predic- 
tions for Af, so we cannot make a detailed co mpar ison 
with our data. But we should point out that 0, 0] 
also obtained much less swelling with / than predicted 
in Re&.jjlUlLHS. 



IV. DISCUSSION 

We have demonstrated that chains growth methods 
with resampling, and the PERM algorithm in particu- 



lar, are able to produce very precise Monte Carlo data 
for star polymers with many arms. Using the Domb- 
Joyce model on the simple cubic lattice, we combined 
this with absence of leading corrections to scaling and 
with the possibility to connect arbitrarily many arms to 
a pointlike core. This allowed us to test conjectured scal- 
ing laws for the entropic critical exponents 7/ and for 
the /-dependent swelling of single arms. In principle we 
could have measured during these simulations also other 
observables like monomer densities, star shapes, radii of 
gyration, etc. 

Our most interesting results are for the exponents 7/. 
All previous simulations were compatible with the predic- 
tions from the heuristic Daoud-Cotton model, but they 
were not very precise. There are also no good experimen- 
tal results for these exponents, although they are funda- 
mental for the entropy (and thus also for the free energy) 
of star polymers in good solvents. Our results show that 
these predictions are qualitatively correct (7/ is negative 
and diverges as — f a , but the exponent a clearly disagrees 
with the prediction. 

We also disagree with the prediction of the Daoud- 
Cotton model for the sizes of star polymers, and indeed 
the disagreement for the end-to-center distances is larger 
than for 7/. They increase with / much slower than pre- 
dicted. But this finding is not entirely new, it had been 
observed previously in Monte Carlo simulations [2(1 |4jj . 
Our data are compatible with these, but more precise 
and extending to larger values of /. Disagreement with 
the Daoud-Cotton prediction for star polymer sizes was 
also found in some experiments l4ll. b ut there are also 
repeated claims in the literature |22j,[42j that experiments 
are compatible with it. We have no good explanation for 
the latter, except that the interpretation of experiments 
for diluted solutions might be less easy than anticipated. 

With slight modifications of the algorithm one can also 
study related problems like stars center-absorbed to sur- 
faces stars confined between two planar walls |44| . 
heterostars uy [40|. interactions between two star poly- 
mers [ill . I38ll45l . l46| . or star polymer - colloid interactions 
|47) . We expect that PERM will be more efficient than 
previous algorithms (not the least because it gives imme- 
diately precise entropy estimates) , in particular if applied 
to lattice models. PERM can also be applied off-lattice 
[3ll |4S[ , but its advantage is in general less pronounced 
there. We hope to present simulations for some of these 
problems in the near future. 
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